home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / moment.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  126 lines
  1. ;$Id: moment.pro,v 1.12 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       MOMENT
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the mean, variance, skewness and kurtosis
  11. ;       of an N-element vector. If the vector contains N identical elements, 
  12. ;       the mean and variance are computed; the skewness and kurtosis are 
  13. ;       not defined and are returned as IEEE NaN (Not a Number). Optionally, 
  14. ;       the mean absolute deviation and standard deviation may also be 
  15. ;       computed.
  16. ;
  17. ; CATEGORY:
  18. ;       Statistics.
  19. ;
  20. ; CALLING SEQUENCE:
  21. ;       Result = Moment(X)
  22. ;
  23. ; INPUTS:
  24. ;       X:      An N-element vector of type integer, float or double.
  25. ;
  26. ; KEYWORD PARAMETERS:
  27. ;       DOUBLE: If set to a non-zero value, computations are done in
  28. ;               double precision arithmetic.
  29. ;
  30. ;       MDEV:   Use this keyword to specify a named variable which returns
  31. ;               the mean absolute deviation of X.
  32. ;
  33. ;       SDEV:   Use this keyword to specify a named variable which returns
  34. ;               the standard deviation of X.
  35. ;
  36. ; EXAMPLE:
  37. ;       Define the N-element vector of sample data.
  38. ;         x = [65, 63, 67, 64, 68, 62, 70, 66, 68, 67, 69, 71, 66, 65, 70]
  39. ;       Compute the mean, variance, skewness and kurtosis.
  40. ;         result = moment(x)
  41. ;       The result should be the 4-element vector: 
  42. ;       [66.7333, 7.06667, -0.0942851, -1.18258]
  43. ;
  44. ; PROCEDURE:
  45. ;       MOMENT computes the first four "moments" about the mean of an N-element
  46. ;       vector of sample data. The computational formulas are given in the IDL 
  47. ;       Reference Guide. 
  48. ;
  49. ; REFERENCE:
  50. ;       APPLIED STATISTICS (third edition)
  51. ;       J. Neter, W. Wasserman, G.A. Whitmore
  52. ;       ISBN 0-205-10328-6
  53. ;
  54. ; MODIFICATION HISTORY:
  55. ;       Written by:  GGS, RSI, August 1994
  56. ;       Modified:    GGS, RSI, September 1995
  57. ;                    Added DOUBLE keyword. 
  58. ;                    Added checking for N identical elements. 
  59. ;                    Added support for IEEE NaN (Not a Number).
  60. ;                    Modified variance formula.
  61. ;       Modified:    GGS, RSI, April 1996
  62. ;                    Modified keyword checking and use of double precision. 
  63. ;-
  64. FUNCTION DATOTAL, Arg, Double = Double
  65.  
  66.   Type = SIZE(Arg)
  67.  
  68.   ;If DATOTAL is called without the Double keyword let the type
  69.   ;of input argument drive the precision of TOTAL.
  70.   if N_ELEMENTS(Double) eq 0 then $
  71.     Double = (Type[Type[0]+1] eq 5) or (Type[Type[0]+1] eq 9)
  72.     
  73.   ArgSum = TOTAL(Arg, Double = Double)
  74.  
  75.   if Type[Type[0]+1] eq 5 and Double eq 0 then $
  76.     RETURN, FLOAT(ArgSum) else $
  77.   if Type[Type[0]+1] eq 9 and Double eq 0 then $
  78.     RETURN, COMPLEX(ArgSum) else $
  79.     RETURN, ArgSum
  80.  
  81. END
  82.  
  83. FUNCTION Moment, X, Double = Double, Mdev = Mdev, Sdev = Sdev
  84.  
  85.   ON_ERROR, 2
  86.  
  87.   TypeX = SIZE(X)
  88.  
  89.   ;Check length.
  90.   if TypeX[TypeX[0]+2] lt 2 then $
  91.     MESSAGE, "X array must contain 2 or more elements."
  92.  
  93.   ;If the DOUBLE keyword is not set then the internal precision and
  94.   ;result are identical to the type of input.
  95.   if N_ELEMENTS(Double) eq 0 then $
  96.     Double = (TypeX[TypeX[0]+1] eq 5 or TypeX[TypeX[0]+1] eq 9)
  97.  
  98.   nX = TypeX[TypeX[0]+2]
  99.   Mean = DATOTAL(X, Double = Double) / nX
  100.   Resid = X - Mean
  101.  
  102.   ;Mean absolute deviation (returned through the Mdev keyword).
  103.   Mdev = DATOTAL(ABS(Resid), Double = Double) / nX
  104.  
  105.   Var = DATOTAL(Resid^2, Double = Double) / (nX-1.0)
  106.     ;Numerically-stable "two-pass" formula.
  107.     ;r2 = DATOTAL(Resid^2, Double = Double)
  108.     ;Var1 = r2 / (nX-1.0)
  109.     ;Var2 = (r2 - (DATOTAL(Resid, Double = Double)^2)/nX)/(nX-1.0)
  110.     ;Var =  (Var1 + Var2)/2.0
  111.  
  112.   ;Standard deviation (returned through the Sdev keyword).
  113.   Sdev = SQRT(Var)
  114.  
  115.   if Sdev ne 0 then begin    ;Skew & kurtosis defined?
  116.     Skew = DATOTAL(Resid^3, Double = Double) / (nX * Sdev ^ 3)
  117.     Kurt = DATOTAL(Resid^4, Double = Double) / (nX * Sdev ^ 4) - 3.0
  118.       ;The "-3" term makes the kurtosis value zero for normal distributions.
  119.       ;Positive values of the kurtosis (lepto-kurtic) indicate pointed or
  120.       ;peaked distributions; Negative values (platy-kurtic) indicate flatt-
  121.       ;ened or non-peaked distributions.
  122.     RETURN, [Mean, Var, Skew, Kurt]
  123.   endif else $        ;All elements equal. Return NaN for skew & kurtosis
  124.     RETURN, [Mean, Var, !VALUES.F_NAN, !VALUES.F_NAN]
  125. end
  126.